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Abstract 

Similar evolutionary variational inequalities appear as convenient for- 
mulations for continuous models for sandpile growth, magnetization of 
type-II superconductors, and evolution of some other dissipative systems 
characterized by the multiplicity of metastable states, long-range inter- 
actions, avalanches, and hysteresis. The origin of this similarity is that 
these are quasistationary models of equilibrium in which the multiplicity 
of metastable states is a consequence of a unilateral condition of equi- 
librium (critical-state constraint). Existing variational formulations for 
critical-state models of sandpiles and superconductors are convenient for 
modelling only the "primary" variables (evolving pile shape and magnetic 
field, respectively). The conjugate variables (the surface sand flux and 
the electric field) are also of interest in various applications. Here we de- 
rive dual variational formulations, similar to mixed variational inequalities 
in plasticity, for the sandpile and superconductor models. These formu- 
lations are used in numerical simulations and allow us to approximate 
simultaneously both the primary and dual variables. 
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Sandpiles and type-II superconductors are examples of spatially extended 
open dissipative systems which have infinitely many metastable states but, 
driven by the external forces, tend to organize themselves into a marginally 
stable "critical state" and are then able to demonstrate almost instantaneous 
long-range interactions. The evolution of such systems is often accompanied by 
sudden collapses, like sandpile avalanches, and hysteresis. Although these are 
dissipative systems of a different nature, their continuous models are equivalent 
to similar variational (or quasivariational) inequalities (see £Q and the references 
therein). The origin of this similarity is that these models are quasistationary 
models of equilibrium and the multiplicity of metastable states is a consequence 
of a unilateral equilibrium condition. The rate with which such a system adjusts 
itself to the changing external conditions is determined implicitly and appears 
in the model as a Lagrange multiplier. Typically, the multiplier is eliminated 
in transition to a variational formulation written in terms of a "primary" vari- 
able (surface of a sandpile, magnetic field in a superconductor, stress tensor in 
elastoplasticity, etc.) In many situations, however, the Lagrange multiplier or, 
equivalcntly, a "dual" variable (sand flux upon the pile surface, electric field, 
and strain tensor, respectively) also has to be found. We present variational 
formulations, for both the sandpile and the superconductivity problem, written 
for the dual variables. Only the simplest version of each problem is considered. 

Sandpiles - Let sand be poured out onto a rigid support surface, y = ho(x), 
given in a domain fl C K 2 with boundary dfl. If the support boundary is open, 
a model for pile surface evolution can be written as 

d t h + V-q = w 1 h\ t =o = h , h\ gn =h \ gn , 
where h is the pile surface, w > is the given source density, q is the horizontal 
projection of the flux of sand pouring down the pile surface. If the support has 
no slopes steeper than the sand angle of repose, |V/io| < 7 = tana, the simplest 
constitutive relations for this model read: (i) the flux q is directed towards the 
steepest decent of the surface, (ii) the surface slope cannot exceed the critical 
angle a, and (iii) the flux is zero upon subcritical slopes. Equivalcntly, one 
can write q = —m\Jh and show that mix, t) > is the Lagrange multiplier 
related to the constraint |V/i| < 7. The model can be rewritten as a variational 
inequality for h, 

h{., t) e K : (d t h - w, <p - h) > V<p G K, h\ t=0 = ho, (1) 
where K = {tp G H 1 ^) : \V(p\ < 7 a.e., f\ 9 Q = h \ dn } and (•, •) is the stan- 
dard L 2 (fl) inner product. Simple analytical solutions of this inequality describe 
piles generated on the support ho = 0. For the point source w = 6(x—Xo), a con- 
ical pile with critical slopes grows until its base touches the domain boundary. 
Then there appears a runway connecting the cone apex with the boundary and 
the pile growth stops: all additional sand just follows the runway and leaves the 
system. On the other hand, if w > everywhere in f2, the final stationary shape 
of the pile is different: h(x) = jdist(x,dil). For w > 0, the general stationary 
solution and an integral representation formula for the corresponding Lagrange 
multiplier to, determining the surface sand flux q, have also been obtained re- 
cently (see and the references therein). Note, however, that it is not easy 
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to compute the Lagrange multiplier using this formula. In the non-stationary 
case, determining the surface flux q remains difficult even if the unique solution 
/i to is found. To compute both these variables simultaneously, we derive a 
dual variational formulation of the evolutionary problem. 

Let {h,q} satisfy the model relations (i)-(iii). Then, for any test field i/>, 
Vh- (V - q) > -\Vh\\il>\ - Wh-q= -|V/i||V| + 7|q| > -7M + 7|<?l- Hence, 
(Vh,ip — q) > 4>{q) — 4>{'4')7 where <f)(q) = 7 J n \q\. Since (V/i,i/> — q) = 
<f gn ho{if> n — q n } — (h, V • {i/> — q}), where ■»/>„ is the normal component of ip 
on dfl, we have </>(i/>) - 4>{q) - (h, V ■ {t/> - q}) + § dn h {ip n - q n } > 0. Noting 
that h = ho + f*wdt — S7- { f Q q dt} we finally obtain 

q(.,t)eV: {^-{l 9*},V-{^-g}^+^-g) + 0(V)-0(g)>O (2) 

for any ip £ V. Here T{u) = § dn hoU n — (Iiq + J* wdt,V ■ u) and we define 
V = {ip € [M(n)] 2 : V ■ ip e L 2 (n)}, where M(0.) is the Banach space of 
bounded Radon measures. 

To approximate (2) numerically, we smoothed the non-diffcrentiablc func- 
tional <p by introducing |q| e = (|q| 2 + e 2 ) 1 / 2 , discretized the regularized equal- 
ity problem in time, employed Raviart-Thomas finite elements of lowest order 
with vertex sampling on the nonlinear term, and solved the resulting nonlin- 
ear algebraic system at each time level iteratively using a form of successive 
over-relaxation (see Fig. ^for an example of a numerical simulation). 




Figure 1: Point source above square support with ho = 0. The cone grows until its base 
touches the support boundary and the runway appears. Shown: 1) sand flux |<j| before and 
after this moment; 2) the final pile shape computed as h = ho + Jg w dt — V ■ {/,* qdt}. 

Superconductors - Phcnomenologically, the magnetic field penetration into 
type-II superconductors can be understood as a nonlinear eddy current problem. 
Let a long cylindrical superconductor with a simply connected cross-section D, 
be placed into a non-stationary uniform external magnetic field h e (t) parallel to 
the cylindrical generators. According to Faraday's law, time variations of this 
field induce in a conductor an electric field e leading to a current j parallel to 
the cross-section plane; this current induces a magnetic field h(x,t) parallel to 
h e . Omitting the displacement current in Maxwell's equations and scaling the 
magnetic permeability to be unity, we obtain the following model, 
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d t (h + h e ) + curie = 0, curl/i=j, h\ t= o = h (x), h\gn = 0, 
where curlw = d Xl U2 — d X2 u\ and curl it = (d X2 u, — d Xl u). Instead of the 
usual Ohm law, a multivalued current-voltage relation (the Bean model) is often 
employed for type-II superconductors. It is postulated that (i) the electric field 
e and the current density j have the same direction, (ii) the current density 
cannot exceed some critical value, j c , and (iii) if the current is subcritical, 
the electric field is zero. One can write e = pj and show that the effective 
resistivity, p(x, t) > 0, is the Lagrange multiplier related to current density 
constraint |j| < j c . Using conditions (i)-(iii) we can eliminate the electric field 
from the model. This yields the variational inequality, 

h(.,t) G K : (d t {h + h e },(p-h)>0 Vtp e K, h\ t=0 = h , 
where K = {tp G Hq(Q,) : |V<p| < 7 a.e.}. This inequality for h can be approx- 
imated numerically and often even solved analytically. However, computing the 
electric field e may be difficult ■ As for the sandpile model, a dual variational 
formulation can be derived to find both variables simultaneously: 

e{.,t)eW: [cm\ edij ,curl{V» - e}j+F(ip-e)+<f>(ip)-(l>(e) > (3) 

for any xp eW, where W = {ip G [M(fl)} 2 : cwclip G L 2 (fl)}, <j>(u) = J Q j c \u\, 
and J-(u) — (h e (t) — h e (0) — ho, curlw). The primary variable, h, is found as 

h = h + h e (0) - h e (t) - curl{ f* edt}. 
The simple transformation R : e = (ei, e'x) — » (e2,— ei) maps W to V and 
enables us to use the same Raviart- Thomas finite element as in the previous 
case. To model the magnetization of a superconductor with a multiply connected 
cross-section (Fig. |2Jl we "filled" the hole and set j c = there. 



e |e| Current streamlines 




Figure 2: Cylindrical superconductor in a growing field (square cross-section with a hole, 
zero initial state). Shown at the same moment in time: the electric field e, level contours of 
|e|, and the current streamlines plotted as levels of h. Note the "runway" (red region in the 
|e| contour plot) through which the magnetic field penetrates the hole and where the electric 
field is the strongest. 
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